##

l = c(0.1,0.25,0.5,1)
b = c(0, 0.2, 1, 5)
prior.mean=0
prior.w=1


fd1 = function(x) l[1]/(1+exp(-b[1]*x)) #attention functions
fd2 = function(x) l[1]/(1+exp(-b[2]*x))
fd3 = function(x) l[1]/(1+exp(-b[3]*x))
fd4 = function(x) l[1]/(1+exp(-b[4]*x))
fd5 = function(x) l[2]/(1+exp(-b[1]*x))
fd6 = function(x) l[2]/(1+exp(-b[2]*x))
fd7 = function(x) l[2]/(1+exp(-b[3]*x))
fd8 = function(x) l[2]/(1+exp(-b[4]*x))
fd9 = function(x) l[3]/(1+exp(-b[1]*x))
fd10 = function(x) l[3]/(1+exp(-b[2]*x))
fd11 = function(x) l[3]/(1+exp(-b[3]*x))
fd12 = function(x) l[3]/(1+exp(-b[4]*x))
fd13 = function(x) l[4]/(1+exp(-b[1]*x))
fd14 = function(x) l[4]/(1+exp(-b[2]*x))
fd15 = function(x) l[4]/(1+exp(-b[3]*x))
fd16 = function(x) l[4]/(1+exp(-b[4]*x))


fr1 = function(x) l[1]/(1+exp(b[1]*x))
fr2 = function(x) l[1]/(1+exp(b[2]*x))
fr3 = function(x) l[1]/(1+exp(b[3]*x))
fr4 = function(x) l[1]/(1+exp(b[4]*x))
fr5 = function(x) l[2]/(1+exp(b[1]*x))
fr6 = function(x) l[2]/(1+exp(b[2]*x))
fr7 = function(x) l[2]/(1+exp(b[3]*x))
fr8 = function(x) l[2]/(1+exp(b[4]*x))
fr9 = function(x) l[3]/(1+exp(b[1]*x))
fr10 = function(x) l[3]/(1+exp(b[2]*x))
fr11 = function(x) l[3]/(1+exp(b[3]*x))
fr12 = function(x) l[3]/(1+exp(b[4]*x))
fr13 = function(x) l[4]/(1+exp(b[1]*x))
fr14 = function(x) l[4]/(1+exp(b[2]*x))
fr15 = function(x) l[4]/(1+exp(b[3]*x))
fr16 = function(x) l[4]/(1+exp(b[4]*x))



d1 <- fd1(seq(from=-3, to = 3, by = 0.01)) #plotting functions
d2 <- fd2(seq(from=-3, to = 3, by = 0.01))
d3 <- fd3(seq(from=-3, to = 3, by = 0.01))
d4 <- fd4(seq(from=-3, to = 3, by = 0.01))
d5 <- fd5(seq(from=-3, to = 3, by = 0.01))
d6 <- fd6(seq(from=-3, to = 3, by = 0.01))
d7 <- fd7(seq(from=-3, to = 3, by = 0.01))
d8 <- fd8(seq(from=-3, to = 3, by = 0.01))
d9 <- fd9(seq(from=-3, to = 3, by = 0.01))
d10 <- fd10(seq(from=-3, to = 3, by = 0.01))
d11 <- fd11(seq(from=-3, to = 3, by = 0.01))
d12 <- fd12(seq(from=-3, to = 3, by = 0.01))
d13 <- fd13(seq(from=-3, to = 3, by = 0.01))
d14 <- fd14(seq(from=-3, to = 3, by = 0.01))
d15 <- fd15(seq(from=-3, to = 3, by = 0.01))
d16 <- fd16(seq(from=-3, to = 3, by = 0.01))


r1 <- fr1(seq(from=-3, to = 3, by = 0.01))
r2 <- fr2(seq(from=-3, to = 3, by = 0.01))
r3 <- fr3(seq(from=-3, to = 3, by = 0.01))
r4 <- fr4(seq(from=-3, to = 3, by = 0.01))
r5 <- fr5(seq(from=-3, to = 3, by = 0.01))
r6 <- fr6(seq(from=-3, to = 3, by = 0.01))
r7 <- fr7(seq(from=-3, to = 3, by = 0.01))
r8 <- fr8(seq(from=-3, to = 3, by = 0.01))
r9 <- fr9(seq(from=-3, to = 3, by = 0.01))
r10 <- fr10(seq(from=-3, to = 3, by = 0.01))
r11 <- fr11(seq(from=-3, to = 3, by = 0.01))
r12 <- fr12(seq(from=-3, to = 3, by = 0.01))
r13 <- fr13(seq(from=-3, to = 3, by = 0.01))
r14 <- fr14(seq(from=-3, to = 3, by = 0.01))
r15 <- fr15(seq(from=-3, to = 3, by = 0.01))
r16 <- fr16(seq(from=-3, to = 3, by = 0.01))


                        
dem <- data.frame(     #storing function plots in a dataframe
  x = seq(-3,3,by=0.01),
  y = c(d1,d2,d3,d4,d5,d6,d7,d8,d9,d10,d11,d12,d13,d14,d15,d16),
  k = factor(c(rep(c(b,b,b,b), 
                   times = c(601,601,601,601,601,601,601,601,601,601,601,601,601,601,601,601)))),
  l = rep(c(l), times = c(2404,2404,2404,2404))
)

dem$b.level=factor(dem$k, levels=c(0,0.2,1,5), labels=c("K=0","K=0.2","K=1","K=5"))
dem$l.level=factor(dem$l, levels=c(0.1,0.25,0.5,1), labels=c("L=0.1","L=0.25","L=0.5","L=1"))




pdem<- ggplot(data = dem, aes(x, y)) +  #creating a figure
  geom_line() + 
  facet_grid(b.level~l.level)+
  theme_bw()+
  theme(strip.background =element_rect(fill="white"))+
  theme(panel.grid.major=element_blank(),
        panel.grid.minor=element_blank())+
  ylab("Attention Rate (%)")+
  xlab("Value of New Information")+
  scale_y_continuous(breaks=c(0,.5, 1),labels=scales::percent)+
  ggtitle("Panel A: Inpartisan's Attention Function")+
  theme(text=element_text(family="Times", size=8))+
  theme(plot.title = element_text(family = "Times", size = (10)))+
  theme(axis.title = element_text(family = "Times", size = (8)))

pdem


rep <- data.frame(
  x = seq(-3,3,by=0.01),
  y = c(r1,r2,r3,r4,r5,r6,r7,r8,r9,r10,r11,r12,r13,r14,r15,r16),
  k = factor(c(rep(c(b,b,b,b), 
                   times = c(601,601,601,601,601,601,601,601,601,601,601,601,601,601,601,601)))),
  l = rep(c(l), times = c(2404,2404,2404,2404))
)
rep$b.level=factor(rep$k, levels=c(0,0.2,1,5), labels=c("K=0","K=0.2","K=1","K=5"))
rep$l.level=factor(rep$l, levels=c(0.1,0.25,0.5,1), labels=c("L=0.1","L=0.25","L=0.5","L=1"))



prep<- ggplot(data = rep, aes(x, y)) + 
  geom_line() + 
  facet_grid(b.level~l.level)+
  theme_bw()+
  theme(strip.background =element_rect(fill="white"))+
  theme(panel.grid.major=element_blank(),
        panel.grid.minor=element_blank())+
  ylab("Attention Rate (%)")+
  xlab("Value of New Information")+
  scale_y_continuous(breaks=c(0,.5, 1),labels=scales::percent)+
  ggtitle("Panel B: Outpartisan's Attention Function")+
  theme(text=element_text(family="Times", size=8))+
  theme(plot.title = element_text(family = "Times", size = (10)))+
  theme(axis.title = element_text(family = "Times", size = (8)))

prep




set.seed (218)
m1=matrix(rep(rnorm(800,mean=0,sd=1)),ncol=100,nrow=8) #simulating "average" info flows
m2=matrix(rep(rnorm(800,mean=-0.5,sd=1)),ncol=100,nrow=8) #simulating "neagtive" info flows

         
d1.post1 = matrix(rep(NA),ncol=100,nrow=9) #posterior belief matrix for each attention function
d1.post2 = matrix(rep(NA),ncol=100,nrow=9)
d1.post3 = matrix(rep(NA),ncol=100,nrow=9)
d1.post4 = matrix(rep(NA),ncol=100,nrow=9)
d1.post5 = matrix(rep(NA),ncol=100,nrow=9)
d1.post6 = matrix(rep(NA),ncol=100,nrow=9)
d1.post7 = matrix(rep(NA),ncol=100,nrow=9)
d1.post8 = matrix(rep(NA),ncol=100,nrow=9)
d1.post9 = matrix(rep(NA),ncol=100,nrow=9)
d1.post10 = matrix(rep(NA),ncol=100,nrow=9)
d1.post11 = matrix(rep(NA),ncol=100,nrow=9)
d1.post12 = matrix(rep(NA),ncol=100,nrow=9)
d1.post13 = matrix(rep(NA),ncol=100,nrow=9)
d1.post14 = matrix(rep(NA),ncol=100,nrow=9)
d1.post15 = matrix(rep(NA),ncol=100,nrow=9)
d1.post16 = matrix(rep(NA),ncol=100,nrow=9)

for (i in 1:100){
d1.post1[,i]=c(cumsum((c(prior.mean,m1[,i]))*c(prior.w,fd1(m1[,i])))/cumsum(c(prior.w,fd1(m1[,i]))))   ###calculating posterior (weighted averages)
d1.post2[,i]=c(cumsum((c(prior.mean,m1[,i]))*c(prior.w,fd2(m1[,i])))/cumsum(c(prior.w,fd2(m1[,i]))))
d1.post3[,i]=c(cumsum((c(prior.mean,m1[,i]))*c(prior.w,fd3(m1[,i])))/cumsum(c(prior.w,fd3(m1[,i]))))
d1.post4[,i]=c(cumsum((c(prior.mean,m1[,i]))*c(prior.w,fd4(m1[,i])))/cumsum(c(prior.w,fd4(m1[,i]))))
d1.post5[,i]=c(cumsum((c(prior.mean,m1[,i]))*c(prior.w,fd5(m1[,i])))/cumsum(c(prior.w,fd5(m1[,i]))))
d1.post6[,i]=c(cumsum((c(prior.mean,m1[,i]))*c(prior.w,fd6(m1[,i])))/cumsum(c(prior.w,fd6(m1[,i]))))
d1.post7[,i]=c(cumsum((c(prior.mean,m1[,i]))*c(prior.w,fd7(m1[,i])))/cumsum(c(prior.w,fd7(m1[,i]))))
d1.post8[,i]=c(cumsum((c(prior.mean,m1[,i]))*c(prior.w,fd8(m1[,i])))/cumsum(c(prior.w,fd8(m1[,i]))))
d1.post9[,i]=c(cumsum((c(prior.mean,m1[,i]))*c(prior.w,fd9(m1[,i])))/cumsum(c(prior.w,fd9(m1[,i]))))
d1.post10[,i]=c(cumsum((c(prior.mean,m1[,i]))*c(prior.w,fd10(m1[,i])))/cumsum(c(prior.w,fd10(m1[,i]))))
d1.post11[,i]=c(cumsum((c(prior.mean,m1[,i]))*c(prior.w,fd11(m1[,i])))/cumsum(c(prior.w,fd11(m1[,i]))))
d1.post12[,i]=c(cumsum((c(prior.mean,m1[,i]))*c(prior.w,fd12(m1[,i])))/cumsum(c(prior.w,fd12(m1[,i]))))
d1.post13[,i]=c(cumsum((c(prior.mean,m1[,i]))*c(prior.w,fd13(m1[,i])))/cumsum(c(prior.w,fd13(m1[,i]))))
d1.post14[,i]=c(cumsum((c(prior.mean,m1[,i]))*c(prior.w,fd14(m1[,i])))/cumsum(c(prior.w,fd14(m1[,i]))))
d1.post15[,i]=c(cumsum((c(prior.mean,m1[,i]))*c(prior.w,fd15(m1[,i])))/cumsum(c(prior.w,fd15(m1[,i]))))
d1.post16[,i]=c(cumsum((c(prior.mean,m1[,i]))*c(prior.w,fd16(m1[,i])))/cumsum(c(prior.w,fd16(m1[,i]))))
}

d1.post= matrix(rep(NA),ncol=16,nrow=900) # reformatting posterior matrices into vectors
d1.post[,1]=as.vector(d1.post1)
d1.post[,2]=as.vector(d1.post2)
d1.post[,3]=as.vector(d1.post3)
d1.post[,4]=as.vector(d1.post4)
d1.post[,5]=as.vector(d1.post5)
d1.post[,6]=as.vector(d1.post6)
d1.post[,7]=as.vector(d1.post7)
d1.post[,8]=as.vector(d1.post8)
d1.post[,9]=as.vector(d1.post9)
d1.post[,10]=as.vector(d1.post10)
d1.post[,11]=as.vector(d1.post11)
d1.post[,12]=as.vector(d1.post12)
d1.post[,13]=as.vector(d1.post13)
d1.post[,14]=as.vector(d1.post14)
d1.post[,15]=as.vector(d1.post15)
d1.post[,16]=as.vector(d1.post16)


d1.result=data.frame(   #data frame storing posterior estimates, time, parameter lists for ggplot setup 
  post=as.vector(d1.post),
t=rep(c(0:8),times=1600),
trial=rep(rep(c(1:100),times=c(rep(9,times=100))),times=16),
b.para = rep(c(rep(b[1],times=900),rep(b[2],times=900),rep(b[3],times=900),rep(b[4],times=900)),times=4),
l.para = rep(c(l), times = c(3600,3600,3600,3600))
)


d1.result$b.level=factor(d1.result$b.para, levels=c(0,0.2,1,5), labels=c("K=0","K=0.2","K=1","K=5"))
d1.result$l.level=factor(d1.result$l.para, levels=c(0.1,0.25,0.5,1), labels=c("L=0.1","L=0.25","L=0.5","L=1"))






d2.post1 = matrix(rep(NA),ncol=100,nrow=9) #data frame storing posterior estimates, time, parameter lists for ggplot setup for negative info flows
d2.post2 = matrix(rep(NA),ncol=100,nrow=9)
d2.post3 = matrix(rep(NA),ncol=100,nrow=9)
d2.post4 = matrix(rep(NA),ncol=100,nrow=9)
d2.post5 = matrix(rep(NA),ncol=100,nrow=9)
d2.post6 = matrix(rep(NA),ncol=100,nrow=9)
d2.post7 = matrix(rep(NA),ncol=100,nrow=9)
d2.post8 = matrix(rep(NA),ncol=100,nrow=9)
d2.post9 = matrix(rep(NA),ncol=100,nrow=9)
d2.post10 = matrix(rep(NA),ncol=100,nrow=9)
d2.post11 = matrix(rep(NA),ncol=100,nrow=9)
d2.post12 = matrix(rep(NA),ncol=100,nrow=9)
d2.post13 = matrix(rep(NA),ncol=100,nrow=9)
d2.post14 = matrix(rep(NA),ncol=100,nrow=9)
d2.post15 = matrix(rep(NA),ncol=100,nrow=9)
d2.post16 = matrix(rep(NA),ncol=100,nrow=9)

for (i in 1:100){
  d2.post1[,i]=c(cumsum((c(prior.mean,m2[,i]))*c(prior.w,fd1(m2[,i])))/cumsum(c(prior.w,fd1(m2[,i]))))   ###weighted averages
  d2.post2[,i]=c(cumsum((c(prior.mean,m2[,i]))*c(prior.w,fd2(m2[,i])))/cumsum(c(prior.w,fd2(m2[,i]))))
  d2.post3[,i]=c(cumsum((c(prior.mean,m2[,i]))*c(prior.w,fd3(m2[,i])))/cumsum(c(prior.w,fd3(m2[,i]))))
  d2.post4[,i]=c(cumsum((c(prior.mean,m2[,i]))*c(prior.w,fd4(m2[,i])))/cumsum(c(prior.w,fd4(m2[,i]))))
  d2.post5[,i]=c(cumsum((c(prior.mean,m2[,i]))*c(prior.w,fd5(m2[,i])))/cumsum(c(prior.w,fd5(m2[,i]))))
  d2.post6[,i]=c(cumsum((c(prior.mean,m2[,i]))*c(prior.w,fd6(m2[,i])))/cumsum(c(prior.w,fd6(m2[,i]))))
  d2.post7[,i]=c(cumsum((c(prior.mean,m2[,i]))*c(prior.w,fd7(m2[,i])))/cumsum(c(prior.w,fd7(m2[,i]))))
  d2.post8[,i]=c(cumsum((c(prior.mean,m2[,i]))*c(prior.w,fd8(m2[,i])))/cumsum(c(prior.w,fd8(m2[,i]))))
  d2.post9[,i]=c(cumsum((c(prior.mean,m2[,i]))*c(prior.w,fd9(m2[,i])))/cumsum(c(prior.w,fd9(m2[,i]))))
  d2.post10[,i]=c(cumsum((c(prior.mean,m2[,i]))*c(prior.w,fd10(m2[,i])))/cumsum(c(prior.w,fd10(m2[,i]))))
  d2.post11[,i]=c(cumsum((c(prior.mean,m2[,i]))*c(prior.w,fd11(m2[,i])))/cumsum(c(prior.w,fd11(m2[,i]))))
  d2.post12[,i]=c(cumsum((c(prior.mean,m2[,i]))*c(prior.w,fd12(m2[,i])))/cumsum(c(prior.w,fd12(m2[,i]))))
  d2.post13[,i]=c(cumsum((c(prior.mean,m2[,i]))*c(prior.w,fd13(m2[,i])))/cumsum(c(prior.w,fd13(m2[,i]))))
  d2.post14[,i]=c(cumsum((c(prior.mean,m2[,i]))*c(prior.w,fd14(m2[,i])))/cumsum(c(prior.w,fd14(m2[,i]))))
  d2.post15[,i]=c(cumsum((c(prior.mean,m2[,i]))*c(prior.w,fd15(m2[,i])))/cumsum(c(prior.w,fd15(m2[,i]))))
  d2.post16[,i]=c(cumsum((c(prior.mean,m2[,i]))*c(prior.w,fd16(m2[,i])))/cumsum(c(prior.w,fd16(m2[,i]))))
}

d2.post= matrix(rep(NA),ncol=16,nrow=900)
d2.post[,1]=as.vector(d2.post1)
d2.post[,2]=as.vector(d2.post2)
d2.post[,3]=as.vector(d2.post3)
d2.post[,4]=as.vector(d2.post4)
d2.post[,5]=as.vector(d2.post5)
d2.post[,6]=as.vector(d2.post6)
d2.post[,7]=as.vector(d2.post7)
d2.post[,8]=as.vector(d2.post8)
d2.post[,9]=as.vector(d2.post9)
d2.post[,10]=as.vector(d2.post10)
d2.post[,11]=as.vector(d2.post11)
d2.post[,12]=as.vector(d2.post12)
d2.post[,13]=as.vector(d2.post13)
d2.post[,14]=as.vector(d2.post14)
d2.post[,15]=as.vector(d2.post15)
d2.post[,16]=as.vector(d2.post16)


d2.result=data.frame(
  post=as.vector(d2.post),
  t=rep(c(0:8),times=1600),
  trial=rep(rep(c(1:100),times=c(rep(9,times=100))),times=16),
  b.para = rep(c(rep(b[1],times=900),rep(b[2],times=900),rep(b[3],times=900),rep(b[4],times=900)),times=4),
  l.para = rep(c(l), times = c(3600,3600,3600,3600))
)


d2.result$b.level=factor(d2.result$b.para, levels=c(0,0.2,1,5), labels=c("K=0","K=0.2","K=1","K=5"))
d2.result$l.level=factor(d2.result$l.para, levels=c(0.1,0.25,0.5,1), labels=c("L=0.1","L=0.25","L=0.5","L=1"))






r1.post1 = matrix(rep(NA),ncol=100,nrow=9)
r1.post2 = matrix(rep(NA),ncol=100,nrow=9)
r1.post3 = matrix(rep(NA),ncol=100,nrow=9)
r1.post4 = matrix(rep(NA),ncol=100,nrow=9)
r1.post5 = matrix(rep(NA),ncol=100,nrow=9)
r1.post6 = matrix(rep(NA),ncol=100,nrow=9)
r1.post7 = matrix(rep(NA),ncol=100,nrow=9)
r1.post8 = matrix(rep(NA),ncol=100,nrow=9)
r1.post9 = matrix(rep(NA),ncol=100,nrow=9)
r1.post10 = matrix(rep(NA),ncol=100,nrow=9)
r1.post11 = matrix(rep(NA),ncol=100,nrow=9)
r1.post12 = matrix(rep(NA),ncol=100,nrow=9)
r1.post13 = matrix(rep(NA),ncol=100,nrow=9)
r1.post14 = matrix(rep(NA),ncol=100,nrow=9)
r1.post15 = matrix(rep(NA),ncol=100,nrow=9)
r1.post16 = matrix(rep(NA),ncol=100,nrow=9)

for (i in 1:100){
  r1.post1[,i]=c(cumsum((c(prior.mean,m1[,i]))*c(prior.w,fr1(m1[,i])))/cumsum(c(prior.w,fr1(m1[,i]))))   ###weighted averages
  r1.post2[,i]=c(cumsum((c(prior.mean,m1[,i]))*c(prior.w,fr2(m1[,i])))/cumsum(c(prior.w,fr2(m1[,i]))))
  r1.post3[,i]=c(cumsum((c(prior.mean,m1[,i]))*c(prior.w,fr3(m1[,i])))/cumsum(c(prior.w,fr3(m1[,i]))))
  r1.post4[,i]=c(cumsum((c(prior.mean,m1[,i]))*c(prior.w,fr4(m1[,i])))/cumsum(c(prior.w,fr4(m1[,i]))))
  r1.post5[,i]=c(cumsum((c(prior.mean,m1[,i]))*c(prior.w,fr5(m1[,i])))/cumsum(c(prior.w,fr5(m1[,i]))))
  r1.post6[,i]=c(cumsum((c(prior.mean,m1[,i]))*c(prior.w,fr6(m1[,i])))/cumsum(c(prior.w,fr6(m1[,i]))))
  r1.post7[,i]=c(cumsum((c(prior.mean,m1[,i]))*c(prior.w,fr7(m1[,i])))/cumsum(c(prior.w,fr7(m1[,i]))))
  r1.post8[,i]=c(cumsum((c(prior.mean,m1[,i]))*c(prior.w,fr8(m1[,i])))/cumsum(c(prior.w,fr8(m1[,i]))))
  r1.post9[,i]=c(cumsum((c(prior.mean,m1[,i]))*c(prior.w,fr9(m1[,i])))/cumsum(c(prior.w,fr9(m1[,i]))))
  r1.post10[,i]=c(cumsum((c(prior.mean,m1[,i]))*c(prior.w,fr10(m1[,i])))/cumsum(c(prior.w,fr10(m1[,i]))))
  r1.post11[,i]=c(cumsum((c(prior.mean,m1[,i]))*c(prior.w,fr11(m1[,i])))/cumsum(c(prior.w,fr11(m1[,i]))))
  r1.post12[,i]=c(cumsum((c(prior.mean,m1[,i]))*c(prior.w,fr12(m1[,i])))/cumsum(c(prior.w,fr12(m1[,i]))))
  r1.post13[,i]=c(cumsum((c(prior.mean,m1[,i]))*c(prior.w,fr13(m1[,i])))/cumsum(c(prior.w,fr13(m1[,i]))))
  r1.post14[,i]=c(cumsum((c(prior.mean,m1[,i]))*c(prior.w,fr14(m1[,i])))/cumsum(c(prior.w,fr14(m1[,i]))))
  r1.post15[,i]=c(cumsum((c(prior.mean,m1[,i]))*c(prior.w,fr15(m1[,i])))/cumsum(c(prior.w,fr15(m1[,i]))))
  r1.post16[,i]=c(cumsum((c(prior.mean,m1[,i]))*c(prior.w,fr16(m1[,i])))/cumsum(c(prior.w,fr16(m1[,i]))))
}

r1.post= matrix(rep(NA),ncol=16,nrow=900)
r1.post[,1]=as.vector(r1.post1)
r1.post[,2]=as.vector(r1.post2)
r1.post[,3]=as.vector(r1.post3)
r1.post[,4]=as.vector(r1.post4)
r1.post[,5]=as.vector(r1.post5)
r1.post[,6]=as.vector(r1.post6)
r1.post[,7]=as.vector(r1.post7)
r1.post[,8]=as.vector(r1.post8)
r1.post[,9]=as.vector(r1.post9)
r1.post[,10]=as.vector(r1.post10)
r1.post[,11]=as.vector(r1.post11)
r1.post[,12]=as.vector(r1.post12)
r1.post[,13]=as.vector(r1.post13)
r1.post[,14]=as.vector(r1.post14)
r1.post[,15]=as.vector(r1.post15)
r1.post[,16]=as.vector(r1.post16)


r1.result=data.frame(
  post=as.vector(r1.post),
  t=rep(c(0:8),times=1600),
  trial=rep(rep(c(1:100),times=c(rep(9,times=100))),times=16),
  b.para = rep(c(rep(b[1],times=900),rep(b[2],times=900),rep(b[3],times=900),rep(b[4],times=900)),times=4),
  l.para = rep(c(l), times = c(3600,3600,3600,3600))
)


r1.result$b.level=factor(r1.result$b.para, levels=c(0,0.2,1,5), labels=c("K=0","K=0.2","K=1","K=5"))
r1.result$l.level=factor(r1.result$l.para, levels=c(0.1,0.25,0.5,1), labels=c("L=0.1","L=0.25","L=0.5","L=1"))






r2.post1 = matrix(rep(NA),ncol=100,nrow=9)
r2.post2 = matrix(rep(NA),ncol=100,nrow=9)
r2.post3 = matrix(rep(NA),ncol=100,nrow=9)
r2.post4 = matrix(rep(NA),ncol=100,nrow=9)
r2.post5 = matrix(rep(NA),ncol=100,nrow=9)
r2.post6 = matrix(rep(NA),ncol=100,nrow=9)
r2.post7 = matrix(rep(NA),ncol=100,nrow=9)
r2.post8 = matrix(rep(NA),ncol=100,nrow=9)
r2.post9 = matrix(rep(NA),ncol=100,nrow=9)
r2.post10 = matrix(rep(NA),ncol=100,nrow=9)
r2.post11 = matrix(rep(NA),ncol=100,nrow=9)
r2.post12 = matrix(rep(NA),ncol=100,nrow=9)
r2.post13 = matrix(rep(NA),ncol=100,nrow=9)
r2.post14 = matrix(rep(NA),ncol=100,nrow=9)
r2.post15 = matrix(rep(NA),ncol=100,nrow=9)
r2.post16 = matrix(rep(NA),ncol=100,nrow=9)

for (i in 1:100){
  r2.post1[,i]=c(cumsum((c(prior.mean,m2[,i]))*c(prior.w,fr1(m2[,i])))/cumsum(c(prior.w,fr1(m2[,i]))))   ###weighted averages
  r2.post2[,i]=c(cumsum((c(prior.mean,m2[,i]))*c(prior.w,fr2(m2[,i])))/cumsum(c(prior.w,fr2(m2[,i]))))
  r2.post3[,i]=c(cumsum((c(prior.mean,m2[,i]))*c(prior.w,fr3(m2[,i])))/cumsum(c(prior.w,fr3(m2[,i]))))
  r2.post4[,i]=c(cumsum((c(prior.mean,m2[,i]))*c(prior.w,fr4(m2[,i])))/cumsum(c(prior.w,fr4(m2[,i]))))
  r2.post5[,i]=c(cumsum((c(prior.mean,m2[,i]))*c(prior.w,fr5(m2[,i])))/cumsum(c(prior.w,fr5(m2[,i]))))
  r2.post6[,i]=c(cumsum((c(prior.mean,m2[,i]))*c(prior.w,fr6(m2[,i])))/cumsum(c(prior.w,fr6(m2[,i]))))
  r2.post7[,i]=c(cumsum((c(prior.mean,m2[,i]))*c(prior.w,fr7(m2[,i])))/cumsum(c(prior.w,fr7(m2[,i]))))
  r2.post8[,i]=c(cumsum((c(prior.mean,m2[,i]))*c(prior.w,fr8(m2[,i])))/cumsum(c(prior.w,fr8(m2[,i]))))
  r2.post9[,i]=c(cumsum((c(prior.mean,m2[,i]))*c(prior.w,fr9(m2[,i])))/cumsum(c(prior.w,fr9(m2[,i]))))
  r2.post10[,i]=c(cumsum((c(prior.mean,m2[,i]))*c(prior.w,fr10(m2[,i])))/cumsum(c(prior.w,fr10(m2[,i]))))
  r2.post11[,i]=c(cumsum((c(prior.mean,m2[,i]))*c(prior.w,fr11(m2[,i])))/cumsum(c(prior.w,fr11(m2[,i]))))
  r2.post12[,i]=c(cumsum((c(prior.mean,m2[,i]))*c(prior.w,fr12(m2[,i])))/cumsum(c(prior.w,fr12(m2[,i]))))
  r2.post13[,i]=c(cumsum((c(prior.mean,m2[,i]))*c(prior.w,fr13(m2[,i])))/cumsum(c(prior.w,fr13(m2[,i]))))
  r2.post14[,i]=c(cumsum((c(prior.mean,m2[,i]))*c(prior.w,fr14(m2[,i])))/cumsum(c(prior.w,fr14(m2[,i]))))
  r2.post15[,i]=c(cumsum((c(prior.mean,m2[,i]))*c(prior.w,fr15(m2[,i])))/cumsum(c(prior.w,fr15(m2[,i]))))
  r2.post16[,i]=c(cumsum((c(prior.mean,m2[,i]))*c(prior.w,fr16(m2[,i])))/cumsum(c(prior.w,fr16(m2[,i]))))
}

r2.post= matrix(rep(NA),ncol=16,nrow=900)
r2.post[,1]=as.vector(r2.post1)
r2.post[,2]=as.vector(r2.post2)
r2.post[,3]=as.vector(r2.post3)
r2.post[,4]=as.vector(r2.post4)
r2.post[,5]=as.vector(r2.post5)
r2.post[,6]=as.vector(r2.post6)
r2.post[,7]=as.vector(r2.post7)
r2.post[,8]=as.vector(r2.post8)
r2.post[,9]=as.vector(r2.post9)
r2.post[,10]=as.vector(r2.post10)
r2.post[,11]=as.vector(r2.post11)
r2.post[,12]=as.vector(r2.post12)
r2.post[,13]=as.vector(r2.post13)
r2.post[,14]=as.vector(r2.post14)
r2.post[,15]=as.vector(r2.post15)
r2.post[,16]=as.vector(r2.post16)


r2.result=data.frame(
  post=as.vector(r2.post),
  t=rep(c(0:8),times=1600),
  trial=rep(rep(c(1:100),times=c(rep(9,times=100))),times=16),
  b.para = rep(c(rep(b[1],times=900),rep(b[2],times=900),rep(b[3],times=900),rep(b[4],times=900)),times=4),
  l.para = rep(c(l), times = c(3600,3600,3600,3600))
)


r2.result$b.level=factor(r2.result$b.para, levels=c(0,0.2,1,5), labels=c("K=0","K=0.2","K=1","K=5"))
r2.result$l.level=factor(r2.result$l.para, levels=c(0.1,0.25,0.5,1), labels=c("L=0.1","L=0.25","L=0.5","L=1"))




pc=ggplot(data=d2.result, aes(x=t, y=post)) +
  geom_hline(yintercept = -.5, linetype="dashed", color="black",alpha=1)+
  geom_line(aes(group=trial),alpha=0.05,colour="red3")+
  geom_smooth(method=loess, se=FALSE, colour="red3", size=0.5) +
  facet_grid(b.level~l.level)+
  theme_bw()+
  theme(strip.background =element_rect(fill="white"))+
  theme(panel.grid.major=element_blank(),
        panel.grid.minor=element_blank())+
  ylab("Perceived Competence")+
  ggtitle("Panel C: Inpartisan's Belief Updating")+
  theme(text=element_text(family="Times", size=8))+
  theme(plot.title = element_text(family = "Times", size = (10)))+
  theme(axis.title = element_text(family = "Times", size = (8)))+
  ylim(-1.5, 1.5)+
  xlab("Belief Updating Rounds")
pc



pd=ggplot(data=r2.result, aes(x=t, y=post)) +
  geom_hline(yintercept = -.5, linetype="dashed", color="black",alpha=1)+
  geom_line(aes(group=trial),alpha=0.05,colour="navy")+
  geom_smooth(method=loess, se=FALSE, colour="navy", size=0.5) +
  facet_grid(b.level~l.level)+
  theme_bw()+
  theme(strip.background =element_rect(fill="white"))+
  theme(panel.grid.major=element_blank(),
        panel.grid.minor=element_blank())+
  ylab("Perceived Competence")+
  ggtitle("Panel D: Outpartisan's Belief Updating")+
  theme(text=element_text(family="Times", size=8))+
  theme(plot.title = element_text(family = "Times", size = (10)))+
  theme(axis.title = element_text(family = "Times", size = (8)))+
  ylim(-1.5, 1.5)+
  xlab("Belief Updating Rounds")
pd




pdf("pc.pdf",width=3,height=3,compress=FALSE)
pc
dev.off()
pdf("pd.pdf",width=3,height=3,compress=FALSE)
pd
dev.off()



pdf("pdem.pdf",width=3,height=3,compress=FALSE)
pdem
dev.off()
pdf("prep.pdf",width=3,height=3,compress=FALSE)
prep
dev.off()
